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Abstract 

The Density Matrix Renormalization Group (DMRG) has become a powerful numerical 
method that can be applied to low- dimensional strongly correlated fermionic and bosonic sys- 
tems. It allows for a very precise calculation of static, dynamical and thermodynamical proper- 
ties. Its field of applicability has now extended beyond Condensed Matter, and is successfully 
used in Statistical Mechanics and High Energy Physics as well. In this article, we briefly review 
the main aspects of the method. We also comment on some of the most relevant applications so 
as to give an overview on the scope and possibilities of DMRG and mention the most important 
extensions of the method such as the calculation of dynamical properties, the application to 
classical systems, inclusion of temperature, phonons and disorder, field theory, time-dependent 
properties and the ab initio calculation of electronic states in molecules. 

1 Introduction 

The basics of the Density Matrix Renormalization Group were developed by S. White in 1992[1] 
and since then DMRG has proved to be a very powerful method for low dimensional interacting 
systems. Its remarkable accuracy can be seen for example in the spin-1 Heisenberg chain: for a 
system of hundreds of sites a precision of 10~ 10 for the ground state energy can be achieved. Since 
then it has been applied to a great variety of systems and problems including, among others, spin 
chains and ladders, fermionic and bosonic systems, disordered models, impurities and molecules and 
2D electrons in high magnetic fields. It has also been improved substantially in several directions 
like two (and three) dimensional (2D) classical systems, stochastic models, the presence of phonons, 
quantum chemistry, field theory, the inclusion of temperature and the calculation of dynamical and 
time-dependent properties. Some calculations have also been performed in 2D quantum systems. 
All these topics are treated in detail and in a pedagogical way in the book [2], where the reader 
can find an extensive review on DMRG. In this article we will attempt to cover the different areas 
where it has been applied without entering into details but in a few cases, where we have chosen 
some representative contributions. We suggest the interested reader to look for further information 
in the referenced work. Our aim here is to give the reader a general overview on the subject. 

One of the most important limitations of numerical calculations in finite systems is the great 
amount of states that have to be considered and its exponential growth with system size. Sev- 
eral methods have been introduced in order to reduce the size of the Hilbert space to be able to 
reach larger systems, such as Monte Carlo, renormalization group (RG) and DMRG. Each method 
considers a particular criterion of keeping the relevant information. 

The DMRG was originally developed to overcome the problems that arise in interacting systems 
in ID when standard RG procedures were applied. Consider a block B (a block is a collection of 



sites) where the Hamiltonian Hb and end-operators are defined. These traditional methods consist 
in putting together two or more blocks (e.g. B-B', which we will call the superblock), connected 
using end-operators, in a basis that is a direct product of the basis of each block, forming Hbb>- This 
Hamiltonian is then diagonalizcd, the superblock is replaced by a new effective block B new formed by 
a certain number m of lowest-lying eigenstates of Hbb' and the iteration is continued (see Ref. [3]). 
Although it has been used successfully in certain cases, this procedure, or similar versions of it, has 
been applied to several interacting systems with poor performance. For example, it has been applied 
to the ID Hubbard model keeping m ~ 1000 states. For 16 sites, an error of 5-10% was obtained 
[4]. Other results [5] were also discouraging. A better performance was obtained [6] by adding a 
single site at a time rather than doubling the block size. However, there is one case where a similar 
version of this method applies very well: the Kondo (and Anderson) model. Wilson[7] mapped the 
one-impurity problem onto a one-dimensional lattice with exponentially descreasing hoppings. The 
difference with the method explained above is that in this case, one site (equivalent to an "onion 
shell" ) is added at each step and, due to the exponential decrease of the hopping, very accurate 
results can be obtained. 

Returning to the problem of putting several blocks together, the main source of error comes 
from the election of eigenstates of Hbb 1 as representative states of a superblock. Since Hbb' has 
no connection to the rest of the lattice, its eigenstates may have unwanted features (like nodes) at 
the ends of the block and this can't be improved by increasing the number of states kept. Based on 
this consideration, Noack and White [8] tried including different boundary conditions and boundary 
strengths. This turned out to work well for single particle and Anderson localization problems but, 
however, it did not improve the results significantly for interacting systems. These considerations led 
to the idea of taking a larger superblock that includes the blocks BB', diagonalize the Hamiltonian 
in this large superblock and then somehow project the most favorable states onto BB'. Then BB' 
is replaced by B new . In this way, awkward features in the boundary would vanish and a better 
representation of the states in the infinite system would be achieved. White[l, 3] proposed the 
density matrix as the optimal way of projecting the best states onto part of the system and this will 
be discussed in the next section. Some considerations concerning the effect of boundary conditions 
in the nature of the states kept (and an analogy to the physics of black holes) is given in [9] . The 
justification of using the density matrix is given in detail in Ref. [2]. A very easy and pedagogical 
way of understanding the basic functioning of DMRG is applying it to the calculation of simple 
quantum problems like one particle in a tight binding chain [10, 11]. 

In the following Section we will briefly describe the standard method; in Sect. 3 we will mention 
some of the most important applications; in Sect. 4 we review the most relevant extensions to the 
method and finally in Sect. 5 we concentrate on the way dynamical calculations can be performed 
within DMRG. 

2 The Method 

The DMRG allows for a systematic truncation of the Hilbert space by keeping the most probable 
states describing a wave function (e.g. the ground state) instead of the lowest energy states usually 
kept in previous real space renormalization techniques. 

The basic idea consists in starting from a small system (e.g with N sites) and then gradually 
increase its size (to N + 2, N + 4,...) until the desired length is reached. Let us call the collection 
of N sites the universe and divide it into two parts: the system and the environment (see Fig. 2). 
The Hamiltonian is constructed in the universe and its ground state |V>o > is obtained. This is 
considered as the state of the universe and called the target state. It has components on the system 
and the environment. We want to obtain the most relevant states of the system, i.e., the states 
of the system that have largest weight in |^ )- To obtain this, the environment is considered as a 
statistical bath and the density matrix[12] is used to obtain the desired information on the system. 
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Figure 1: A scheme of the superblock (universe) configuration for the DMRG algorithm[3]. 



So instead of keeping eigenstates of the Hamiltonian in the block (system), we keep eigenstates of 
the density matrix. We will be more explicit below. 

Let's define block [B] as a finite chain with I sites, having an associated Hilbcrt space with, m 
states where operators are defined (in particular the Hamiltonian in this finite chain, Hb and the 
operators at the ends of the block, useful for linking it to other chains or added sites). Except 
for the first iteration, the basis in this block isn't explicitly known due to previous basis rotations 
and reductions. The operators in this basis are matrices and the basis states are characterized by 
quantum numbers (like S z , charge or number of particles, etc). We also define an added block or 
site as [a] having n states. A general iteration of the method is described below: 

i) Define the Hamiltonian Hbb' for the superblock (the universe) formed by putting together 
two blocks [B] and [B'] and two added sites [a] and [a'] in this way: [B a a' B' ] (the primes 
are only to indicate additional blocks, but the primed blocks have the same structure as the non- 
primed ones; this can vary, see the finite-size algorithm below). In general, blocks [B] and [B'] come 
from the previous iteration. The total Hilbert space of this superblock is the direct product of the 
individual spaces corresponding to each block and the added sites. In practice a quantum number of 
the superblock can be fixed (in a spin chain for example one can look at the total S z = subspace) , 
so the total number of states in the superblock is much smaller than (mn) 2 . In some cases, as the 
quantum number of the superblock consists of the sum of the quantum numbers of the individual 
blocks, each block must contain several subspaces (several values of S z for example). Here periodic 
boundary conditions can be attached to the ends and a different block layout should be considered 
(e.g. [B a B' a' ]) to avoid connecting blocks [B] and [B'] which takes longer to converge. The 
boundary conditions are between [a'] and [B]. For closed chains the performance is poorer than for 
open boundary conditions [3, 13]. 

ii) Diagonalize the Hamiltonian Hbb 1 to obtain the ground state | Vo) (target state) using 
Lanczos[14] or Davidson[15] algorithms. Other states could also be kept, such as the first excited 
ones: they are all called target states. 

iii) Construct the density matrix: 

pw = ^2ipo,ijipo,i'j (1) 

3 

on block [B a], where i/)o,ij — (i ® j\i>o) , the states \i) belonging to the Hilbert space of the block [B 
a] and the states \j) to the block [B' a']. The density matrix considers the part [B a] as a system 
and [B' a'], as a statistical bath. The eigenstates of p with the highest eigenvalues correspond 
to the most probable states (or equivalently the states with highest weight) of block [B a] in the 
ground state of the whole superblock. These states are kept up to a certain cutoff, keeping a total 
of m states per block. The density matrix eigenvalues sum up to unity and the truncation error, 
defined as the sum of the density matrix eigenvalues corresponding to discarded eigenvectors, gives 
a qualitative indication of the accuracy of the calculation. 
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iv) With these m states a rectangular matrix O is formed and it is used to change basis and 
reduce all operators defined in [B a]. This block [B a] is then renamed as block [B neMJ ] or simply 
[B] (for example, the Hamiltonian in block [B a], Hsa, is transformed into Hb as Hb = O^HsaO). 

v) A new block [a] is added (one site in our case) and the new superblock [B a a' B'] is formed 
as the direct product of the states of all the blocks. 

vi) This iteration continues until the desired length is achieved. At each step the length is 
N = 21 + 2 (if [a] consists of one site) . 

When more than one target state is used, i.e more than one state is wished to be well described, 
the density matrix is defined as: 



where pi defines the probability of finding the system in the target state \cj>i} (not necessarily eigen- 
statcs of the Hamiltonian). 

The method described above is usually called the infinite- system algorithm since the system size 
increases at each iteration. There is a way to increase precision at each length N called the finite- 
system algorithm. It consists of fixing the lattice size and zipping a couple of times until convergence 
is reached. In this case and for the block configuration [B a a' B' ], N = 1 + 1 + 1 + 1' where I 
and I' are the number of sites in B and B' respectively. In this step the density matrix is used to 
project onto the left I + 1 sites. In order to keep N fixed, in the next block configuration, the right 
block B' should be defined in I - 1 sites such that N={1 + 1) + 1 + 1 + {1- I)'. The operators in 
this smaller block should be kept from previous iterations (in some cases from the iterations for the 
system size with N — 2) [2]. 

The calculation of static properties like correlation functions is easily done by keeping the oper- 
ators in question at each step and performing the corresponding basis change and reduction, in a 
similar manner as done with the Hamiltonian in each block[3]. The energy and measurements are 
calculated in the superblock. In Rcf. [16] an interpretation of the correlation functions of systems at 
criticality is given in terms of wave function entanglement, conjecturing a modification of DMRG 
for these cases that preserves the entanglement. 

A faster convergence of Lanczos or Davidson algorithm is achieved by choosing a good trial 
vector[I7, 18]. An interesting analysis on DMRG accuracy is done in Ref. [19]. Fixed points of the 
DMRG and their relation to matrix product wave functions were studied in [20] and an analytic 
formulation combining the block renormalization group with variational and Fokker-Planck methods 
in [21]. The connection of the method with quantum groups and conformal field theory is treated 
in [22]. There are also interesting connections between the density matrix spectra and integrable 
models [23] via corner transfer matrices. These articles give a deep insight into the essence of the 
DMRG method. 

3 Applications 

Since its development, the number of papers using DMRG has grown enormously and other improve- 
ments to the method have been performed. We would like to mention some applications where this 
method has proved to be useful. Other applications related to further developments of the DMRG 
will be mentioned in Sect. 4. 

A very impressive result with unprecedented accuracy was obtained by White and Huse [24] when 
calculating the spin gap in a S — 1 Hciscnbcrg chain obtaining A = 0.41050J. They also calculated 
very accurate spin correlation functions and excitation energies for one and several magnon states 
and performed a very detailed analysis of the excitations for different momenta. They obtained a spin 
correlation length of 6.03 lattice spacings. Simultaneously S0rcnscn and Affleck [25] also calculated 
the structure factor and spin gap for this system up to length 100 with very high accuracy, comparing 
their results with the nonlinear a model. In a subsequent paper [26] they applied the DMRG to the 
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anisotropic 5 = 1 chain, obtaining the values for the Haldane gap. They also performed a detailed 
study of the S —1/2 end excitations in an open chain. Thermodynamical properties in open S = 1 
chains such as specific heat, electron paramagnetic resonance (EPR) and magnetic susceptibility 
calculated using DMRG gave an excellent fit to experimental data, confirming the existence of free 
spins 1/2 at the boundaries [27]. A related problem, i.e. the effect of non-magnetic impurities 
in spin systems (dimerized, ladders and 2D) was studied in [28, 29]. In addition, the study of 
magnon interactions and magnetization of 5 — 1 chains was done in [30], supersymmetric spin 
chains modelling plateau transitions in the integer quantum Hall effect in [31] and ESR studies in 
these systems was considered in [32]. For larger integer spins there have also been some studies. 
Nishiyama and coworkers [33] calculated the low energy spectrum and correlation functions of the 
5 = 2 antiferromagnetic Heisenberg open chain. They found 5=1 end excitations (in agreement 
with the Valence Bond Theory). Edge excitations for other values of 5 have been studied in Ref. 
[34]. Almost at the same time Schollwock and Jolicoeur[35] calculated the spin gap in the same 
system, up to 350 sites, (A = 0.085J), correlation functions that showed topological order and a 
spin correlation length of 49 lattice spacings. More recent accurate studies of 5 = 2 chains are found 
in [36, 37, 38] and of 5 — 1 chains in staggered magnetic fields [39] including a dctalicd comparison 
to the non-linear sigma model in [40]. In Ref. [41] the dispersion of the single magnon band and 
other properties of the 5 = 2 antiferromagnetic Heisenberg chains were calculated. 

Concerning 5=1/2 systems, DMRG has been crucial for obtaining the logarithmic corrections to 
the 1 /r dependence of the spin-spin correlation functions in the isotropic Heisenberg model [42] . For 
this, very accurate values for the energy and correlation functions were needed. For N = 100 sites 
an error of 10 -5 was achieved keeping m = 150 states per block, comparing with the exact finite-size 
Bethe Ansatz results. For this model it was found that the data for the correlation function has a 
very accurate scaling behaviour and advantage of this was taken to obtain the logarithmic corrections 
in the thermodynamic limit. Other calculations of the spin correlations have been performed for the 
isotropic [43, 44] and anisotropic cases [45]. Luttinger liquid behaviour with magnetic fields have 
been studied in [46], field- induced gaps in [47], anisotropic systems in [48, 49] and the Heisenberg 
model with a weak link in [50] . An analysis of quantum critical points and critical behaviour in spin 
chains by combining DMRG with finite-size scaling was done in [51]. 

Similar calculations have been performed for the 5 = 3/2 Heisenberg chain [52]. In this case a 
stronger logarithmic correction to the spin correlation function was found. For this model there was 
interest in obtaining the central charge c to elucidate whether this model corresponds to the same 
universality class as the 5=1/2 case, where the central charge can be obtained from the finite-size 
scaling of the energy. Although there have been previous attempts [53], these calculations presented 
difficulties since they involved also a term ~ 1/ In 3 N. With the DMRG the value c—1 was clearly 
obtained. 

In Ref. [54], DMRG was applied to an effective spin Hamiltonian obtained from an SU(4) spin- 
orbit critical state in ID. Other applications were done to enlarged symmetry cases with SU(4) 
symmetry in order to study coherence in arrays of quantum dots[55], to obtain the phase diagram 
for ID spin orbital models[56] and dynamical properties in a magnetic field [5 7]. 

Dimerization and frustration have been considered in Refs. [58, 59, 60, 61, 62, 63, 64, 65, 66] 
and alternating spin chains in [67]. 

The case of several coupled spin chains (ladder models) have been investigated in [68, 69, 70, 71, 
72], spin ladders with cyclic four-spin exchanges in [73, 74, 75, 76] and Kagome antiferromagncts 
in [77]. Zigzag spin chains have been considered in [78, 79, 80] and spin chains of coupled triangles 
in [81, 82, 83]. As the DMRG's performance is optimal in open systems, an interesting analysis of 
the boundary effect on correlation functions is done in [13]. Magnetization properties and plateaus 
for quantum spin ladder systems [84, 85, 86] have also been studied. An interesting review on the 
applications to some exact and analytical techniques for quantum magnetism in low dimension, 
including DMRG, is presented in [87]. 

There has been a great amount of applications to fermionic systems such as ID Hubbard and t-J 
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models [88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98], Luttinger liquids with boundaries [99], the Falicov- 
Kimball model [100], the quasiperiodic Aubry- Andre chain[101] and Fibonacci-Hubbard modcls[102]. 
It has also been applied to field theory [9, 103]. The method has been very successful for several band 
Hubbard models[104], Hubbard ladders [105, 106, 107] and t-J ladders[108]. Also several coupled 
chains at different dopings have been considered [109, 110] as well as flux phases in these systems 
[111]. Time reversal symmetry-broken fermionic ladders have been studied in [112] and power laws 
in spinless-fermion ladders in [113]. Long-range Coulomb interactions in the ID electron gas and 
the formation of a Wigner crystal was studied in [114]. Several phases including the Wigner crystal, 
incompressible and compressible liquid states, stripe and pairing phases, have been found using 
DMRG for 2D electrons in high magnetic fields considering different Landau levels[115]. Persistent 
currents in mesoscopic systems have been considered in [116]. 

Quite large quasi-2D systems can be reached, for example in [117] where a 4x20 lattice was 
considered to study ferromagnetism in the infinitc-U Hubbard model; the ground state of a 4-leg t-J 
ladder in [118]; the one and two hole ground state in 9x9 and 10x7 t-J lattices in [119]; a doped 3-leg 
t-J ladder in [120]; the study of striped phases in [121]; domain walls in 19x8 t-J systems in [122]; 
the 2D t-J model in [123] and the magnetic polaron in a 9x9 t-J lattice in [124]. Also big CaVtJDo, 
spin-1/2 lattices reaching 24x11 sites[17] have been studied. There have been some recent attempts 
to implement DMRG in two and higher dimensions [125, 126, 127, 128, 129] but the performance is 
still poorer than in ID. A recent extension using a two-step DMRG algorithm for highly anisotropic 
spin systems has shown promising results [130]. 

Impurity problems have been studied for example in one- [131] and two-impurity [132] Kondo 
systems, in spin chains [133] and in Luttinger Liquids [134]. There have also been applications to 
Kondo and Anderson lattices [135, 136, 137, 138, 139, 140, 141, 142, 143, 144], Kondo lattices with 
localized f 2 configurations [145], the two-channel Kondo lattice on a ladder [146], a t-J chain coupled 
to localized Kondo spins[147] and ferromagnetic Kondo models for manganites [148, 149, 150]. 

4 Other extensions to DMRG 

There have been several extensions to DMRG like the inclusion of symmetries to the method such 
as spin and parity [151, 152, 153]. Total spin conservation and continuous symmetries have been 
treated in [143] and in interaction-round a face Hamiltonians[154], a formulation that can be ap- 
plied to rotational- invariant sytems like S — 1 and 2 chains [37]. A momentum representation of 
this technique [110, 155, 126] that allows for a diagonalization in a fixed momentum subspace has 
been developed as well as applications in dimension higher than one[17, 125, 126, 156] and Bcthc 
latticcs[157]. The inclusion of symmetries is essential to the method since it allows to consider a 
smaller number of states, enhance precision and obtain eigenstates with definite quantum numbers. 
Other recent applications have been in nuclear shell model calculations where a two level pairing 
model has been considered [158] and in the study of ultrasmall superconducting grains, in this case, 
using the particle (hole) states around the Fermi level as the system (environment) block[159]. 

A very interesting and successful application is a recent work in High Energy Physics[160]. Here 
the DMRG is used in an asymptotically free model with bound states, a toy model for quantum 
chromodynamics, namely the two dimensional delta-function potential. For this case an algorithm 
similar to the momentum space DMRG [155] was used where the block and environment consist of 
low and high energy states respectively. The results obtained here are much more accurate than with 
the similarity renormalization group[161] and a generalization to field-theoretical models is proposed 
based on the discreet light-cone quantization in momentum space[162]. Below we briefly mention 
other important extensions, leaving the calculation of dynamical properties for the next Section. 
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4.1 Classical systems 

The DMRG has been very successfully extended to study classical systems. For a detailed description 
we refer the reader to Ref. [163]. Since ID quantum systems are related to 2D classical systems[164], 
it is natural to adapt DMRG to the classical 2D case. This method is based on the renormalization 
group transformation for the transfer matrix T (TMRG). It is a variational method that maximizes 
the partition function using a limited number of degrees of freedom, where the variational state is 
written as a product of local matrices [20]. For 2D classical systems, this algorithm is superior to 
the classical Monte Carlo method in accuracy, speed and in the possibility of treating much larger 
systems. A recent improvement of this method considering periodic boundary conditions is given 
in [165] and a detailed comparison between symmetric and asymmetric targetting is done in [166]. 
TMRG has also been successfully used to renormalize stochastic transfer matrices in a study of 
cellular automatons[167]. The calculation of thermodynamical properties of 3D classical statistical 
systems has been proposed[127] where the eigenstate of the transfer matrix with maximum eigenvalue 
is represented by the product of local tensors optimized using DMRG. 

A further improvement to this method is based on the corner transfer matrix[168], the CTMRG[169, 
170, 171, 172] and can be generalized to any dimension [173]. 

It was first applied to the Ising modcl[163, 174, 175, 176] and also to the Potts modcl[177], 
where very accurate density profiles and critical indices were calculated. Further applications have 
included non-hermitian problems in equilibrium and non-equilibrium physics. In the first case, 
transfer matrices may be non-hermitian and several situations have been considered: a model for the 
Quantum Hall effect [178], the g-symmetric Heisenberg chain related to the conformal series of critical 
models[179] and the anisotropic triangular nearest and next-nearest neighbour Ising models [83]. 
In the second case, the adaptation of the DMRG to non-equilibrium physics like the asymmetric 
exclusion problcm[180] and reaction-diffusion problems [181, 182] has shown to be very successful. 
It has also been applied to stochastic lattice models like in [183] and to the 2D XY model [184]. 

4.2 Finite-temperature DMRG 

The adaptation of the DMRG method for classical systems paved the way for the study of ID 
quantum systems at non zero temperature, by using the Trotter-Suzuki method [185, 186, 187, 188, 
189]. In this case the system is infinite and the finiteness is in the level of the Trotter approximation. 
Standard DMRG usually produces its best results for the ground state energy and less accurate 
results for higher excitations. A different situation occurs here: the lower the temperature, the less 
accurate the results. 

Very nice results have been obtained for the dimerized, S = 1/2, XY model, where the specific 
heat was calculated involving an extremely small basis set [185] (m = 16), the agreement with 
the exact solution being much better in the case where the system has a substantial gap. It has 
also been used to calculate thermodynamical properties of the anisotropic S = 1/2 Heisenberg 
model, with relative errors for the spin susceptibility of less than 10 -3 down to temperatures of 
the order of 0.01J keeping m — 80 states[187]. A complete study of thermodynamical properties 
like magnetization, susceptibility, specific heat and temperature dependent correlation functions 
for the S = 1/2 and 3/2 Heisenberg models was done in [190]. Other applications have been the 
calculation of the temperature dependence of the charge and spin gap in the Kondo insulator [191], 
the calculation of thermodynamical properties of fcrrimagnetic chains[192] and spin ladders[86], the 
study of impurity properties in spin chains[193, 194], frustrated quantum spin chains[195], t-J[196] 
and spin ladders[197] and dimerized frustrated Heisenberg chains[198]. 

An alternative way of incorporating temperature into the DMRG procedure was developed by 
Moukouri and Caron[199]. They considered the standard DMRG taking into account several low- 
lying target states (see Eq. 2) to construct the density matrix, weighted with the Boltzmann factor 
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(fi is the inverse temperature): 

l 3 

With this method they performed reliable calculations of the magnetic susceptibility of quantum spin 
chains with S =1/2 and 3/2, showing excellent agreement with Bethe Ansatz exact results. They 
also calculated low temperature thermodynamical properties of the ID Kondo Lattice Modcl[200] 
and of organic conductors [201]. Zhang et al.[202] applied the same method in the study of a 
magnetic impurity embedded in a quantum spin chain. 

4.3 Phonons, bosons and disorder 

A significant limitation to the DMRG method is that it requires a finite basis and calculations in 
problems with infinite degrees of freedom per site require a large truncation of the basis states [203]. 
However, Jeckelmann and White developed a way of including phonons in DMRG calculations 
by transforming each boson site into several artificial interacting two-state pseudo-sites and then 
applying DMRG to this interacting system [204] (called the "pseudo-site system"). The idea is based 
on the fact that DMRG is much better able to handle several few-states sites than few many-state 
sites [205]. The key idea is to substitute each boson site with 2 N states into N pseudo-sites with 
2 states [206]. They applied this method to the Holstcin model for several hundred sites (keeping 
more than a hundred states per phonon mode) obtaining negligible error. In addition, up to date, 
this method is the most accurate one to determine the ground state energy of the polaron problem 
(Holstein model with a single electron). 

An alternative method (the "Optimal phonon basis" ) [207] is a procedure for generating a con- 
trolled truncation of a large Hilbert space, which allows the use of a very small optimal basis without 
significant loss of accuracy. The system here consists of only one site and the environment has sev- 
eral sites, both having electronic and phononic degrees of freedom. The density matrix is used to 
trace out the degrees of freedom of the environment and extract the most relevant states of the site 
in question. In following steps, more bare phonons are included to the optimal basis obtained in 
this way. This method was successfully applied to study the interactions induced by quantum fluc- 
tuations in quantum strings, as an application to cuprate stripes[208] and the dissipative two-state 
system [209]. A variant of this scheme is the "four block method", as described in [210]. They obtain 
very accurately the Luttinger liquid-CDW insulator transition in the ID Holstein model for spinless 
fermions. 

The method has also been applied to pure bosonic systems such as the disordered bosonic Hub- 
bard modcl[211], where gaps, correlation functions and superfluid density are obtained. The phase 
diagram for the non-disordered Bosc-Hubbard model, showing a reentrance of the superfluid phase 
into the insulating phase was calculated in Ref. [212]. It has also been used to study a chain of 
oscillators with optical phonon spectrum [213] and optical phonons in the quarter-filled Hubbard 
model for organic conductors [2 14]. 

The DMRG has also been generalized to ID random and disordered systems, and applied to 
the random antiferromagnetic and ferromagnetic Heisenberg chains [2 15], including quasiperiodic 
exchange modulation[216] and a detailed study of the Haldane[217] and Griffiths phase[218] in these 
systems. Strongly disordered spin ladders have been considered in [219]. It has also been used in 
disordered Fermi systems such as the spinless modcl[220, 221]. In particular, the transition from the 
Fermi glass to the Mott insulator and the strong enhancement of persistent currents in the transition 
was studied in correlated one-dimensional disordered rings [222]. Disorder-induced crossover effects 
at quantum critical points were studied in [223]. 
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4.4 Molecules and Quantum Chemistry 

There have been several applications to molecules and polymers, such as the Pariser-Parr-Pople 
(PPP) Hamiltonian for a cyclic polyene[224] (where long-range interactions are included), magnetic 
Keplerate molecules [225], molecular Iron rings [226] and polyacenes considering long range inter- 
actions [227]. It has also been applied to conjugated organic systems (polymers), adapting the 
DMRG to take into account the most important symmetries in order to obtain the desired excited 
states[151]. Also conjugated one-dimensional semiconductors [228] have been studied, in which the 
standard approach can be extended to complex ID oligomers where the fundamental repeat is not 
just one or two atoms, but a complex molecular building block. Relatively new fields of application 
are the calculation of dynamical properties in the Rubinstein-Duke model for reptons [229] and 
excitons in dendrimcr molecules [230]. 

Recent attempts to apply DMRG to the ab initio calculation of electronic states in molecules 
have been successful [23 1, 232, 233, 234]. Here, DMRG is applied within the conventional quantum 
chemical framework of a finite basis set with non-orthogonal basis functions centered on each atom. 
After the standard Hartree-Fock (HF) calculation in which a Hamiltonian is produced within the 
orthogonal HF basis, DMRG is used to include correlations beyond HF, where each orbital is treated 
as a "site" in a ID lattice. One important difference with standard DMRG is that, as the interactions 
are long-ranged, several operators must be kept, making the calculation somewhat cumbersome. 
However, very accurate results have been obtained in a check performed in a water molecule (keeping 
up to 25 orbitals and m ~ 200 states per block), obtaining an offset of 0.00024Hartrees with respect to 
the exact ground state energy[235], a better performance than any other approximate method[231]. 

In order to avoid the non-locality introduced in the treatment explained above, White introduced 
the concept of orthlets, local, orthogonal and compact wave functions that allow prior knowledge 
about singularities to be incorporated into the basis and an adequate resolution for the cores[232]. 
The most relevant functions in this basis are chosen via the density matrix. An application based 
on the combination with the momentum version of DMRG is used in [236] to calculate the ground 
state of several molecules. 

5 Dynamical correlation functions 

The DMRG was originally developed to calculate static ground state properties and low-lying ener- 
gies. However, it can also be useful to calculate dynamical response functions. These are of great 
interest in condensed matter physics in connection with experiments such as nuclear magnetic res- 
onance (NMR), neutron scattering, optical absorption, photoemission, etc. We will describe three 
different methods in this Section. 

A recent development for calculating response functions in single impurity systems in the presence 
of a magnetic field was done in [237] by using the DMRG within Wilson's NRG to obtain the Green's 
function. 

An interesting extension of DMRG to tackle time-dependent quantum many-body systems out 
of equilibrium was considered by Cazalilla and Marston[238]. 

5.1 Lanczos and correction vector techniques 

An effective way of extending the basic ideas of this method to the calculation of dynamical quan- 
tities is described in Ref.[239]. It is important to notice here that due to the particular real-space 
construction, it is not possible to fix the momentum as a quantum number. However, we will show 
that by keeping the appropriate target states, a good value of momentum can be obtained. 
We want to calculate the following dynamical correlation function at T = 0: 

C A (t-t') = {^(t)A(t')\^ ), (4) 
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where is the Hermitean conjugate of the operator A, A(t) is the Heisenberg representation of A, 
and \ipo) is the ground state of the system. Its Fourier transform is: 



CUM - E \&n\A\%)\ 2 5(lu - (E n - E )), (5) 

n 

where the summation is taken over all the eigenstates \tp n ) of the Hamiltonian H with energy E n , 
and Eo is the ground state energy. 
Defining the Green's function 

G A (z) = (^\A\z-H)- 1 A\^) 1 (6) 
the correlation function CU(w) can be obtained as 

Ca{u) = -- lim ImGifw + irz + Bo). (7) 
The function Ga can be written in the form of a continued fraction: 

Ga[z) _ MAUM (8) 

z - a "-js- 

z — a-i — g _ 2 

The coefficients a„ and 6„ can be obtained using the following recursion equations [240, 241]: 

\f n+1 )=H\f n )-a n \f n )-b 2 n \f n . 1 ) (9) 

where 

l/o) = A\^ ) 

an = (fn\H\f n )/{f n \f n ), 

(fn\fn)/(fn-l\fn-lh = (10) 



(,2 



For finite systems the Green's function G A (z) has a finite number of poles so only a certain 
number of coefficients a n and b n have to be calculated. The DMRG technique presents a good 
framework to calculate such quantities. With it, the ground state, Hamiltonian and the operator A 
required for the evaluation of CU(w) are obtained. An important requirement is that the reduced 
Hilbert space should also describe with great precision the relevant excited states This is 

achieved by choosing the appropriate target states. For most systems it is enough to consider as 
target states the ground state |^o) and the first few with n = 0, 1... and |/ ) = A\i(j ) as described 
above. In doing so, states in the reduced Hilbert space relevant to the excited states connected to 
the ground state via the operator of interest A are included. The fact that |/o) is an excellent trial 
state, in particular, for the lowest triplet excitations of the two-dimensional antiferromagnet was 
shown in Rcf. [242]. Of course, if the number m of states kept per block is fixed, the more target 
states considered, the less precisely each one of them is described. An optimal number of target 
states and m have to be found for each case. Due to this reduction, the algorithm can be applied up 
to certain lengths, depending on the states involved. For longer chains, the higher energy excitations 
will become inaccurate. Proper sum rules have to be calculated to determine the errors in each case. 

As an application of the method we calculate 

S"(q,w) = ]T |(Vn|^|Vo}| 2 6{w - (E n - E Q )), (11) 

n 

for the ID isotropic Heisenberg model with spin S = 1/2. 
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The spin dynamics of this model has been extensively studied. The lowest excited states in the 
thermodynamic limit are the des Cloiseaux-Pearson triplets [243] , having total spin S T = 1 . The 
dispersion of this spin-wave branch is u) l g = 4f|sin(g)|. Above this lower boundary there exists 
a two-parameter continuum of excited triplet states that have been calculated using the Bcthc 
Ansatz approach [244] with an upper boundary given by u)g — Jn\ sin(g/2)|. It has been shown 
[245], however, that there are excitations above this upper boundary due to higher order scattering 
processes, with a weight that is at least one order of magnitude lower than the spin-wave continuum. 

In Fig. 2 we show the spectrum for q — tt and N — 24 for different values of m, where exact 
results are available for comparison. The delta peaks of Eq. (11) are broadened by a Lorentzian 
for visualizing purposes. As expected, increasing m gives more precise results for the higher excita- 
tions. This spectra has been obtained using the infinite-system method and more precise results are 
expected using the finite-system method, as described later. 




Figure 2: Spectral function for a Heiscnberg chain with TV = 24 and q = it. Full line: exact result 
[246]. The rest are calculated using DMRG with m = 100 (long-dashed line), m = 150 (dashed line) 
and m = 200 (dotted line). 

In Fig. 3 we show the spectrum for two systems lengths and q — tt and q = it/2 keeping m = 200 
states and periodic boundary conditions. For this case it was enough to take 3 target states, i. 
e - | "00 } 7 l/o) = and \fi). Here we have used ~ 40 pairs of coefficients a n and b n , but we 

noticed that if we considered only the first (~ 10) coefficients a n and b n , the spectrum at low energies 
remains essentially unchanged. Minor differences arise at ui/J ~ 2. This is another indication that 
only the first |/„) are relevant for the low energy dynamical properties for finite systems. 

In the inset of Fig. 3 the spectrum for q = 7r/2 and N ~ 28 is shown. For this case we considered 
5 target states i. e. \ipo), |/o) = S*/ 2 \ipo), \f n ) n — 1,3 and m — 200. Here, and for all the cases 
considered, we have verified that the results are very weakly dependent on the weights pi of the 
target states (see Eq.(2)) as long as the appropriate target states are chosen. For lengths where this 
value of q is not defined we took the nearest value. 

Even though we are including states with a given momentum as target states, due to the particular 
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Figure 3: Spectral densities for q = ir, N = 28 (continuous line) and N = 40 (dotted line). Inset: 
Spectral density for q = tt/2 for N = 28 [r] = 0.05). 

real-space construction of the reduced Hilbert space, this translational symmetry is not fulfilled 
and the momentum is not fixed. To check how the reduction on the Hilbert space influences the 
momentum q of the target state |/o) = S*\ipo), we calculated the expectation values (ipo\S*L q , Sg\ipo) 
for all q'. If the momenta of the states were well defined, this value is proportional to S q - q < if q ^ 0. 
fbrg = 0,£ r 3f = 0. 

The momentum distribution for q — ir is shown in Fig. 4 in a semilogarithmic scale where the 
y-axis has been shifted by .003 so as to have well-defined logarithms. We can see here that the 
momentum is better defined, even for much larger systems, but, as expected, more weight on other 
q' values arises for larger N. 

As a check of the approximation we calculated the sum rule 

-i poo p 2tt i 

^ ^ J S"{q,u>) = (M(S z r=a f\^o) = i (12) 

for N = 28, 5 target states and m = 200. We obtain a relative error of 0.86%. 

Recently, important improvements to this method have been published [247] : By considering the 
finite system method in open chains, Kiihner and White obtained a higher precision in dynamical 
responses of spin chains. In order to define a momentum in an open chain and to avoid end effects, 
they introduce a filter function with weight centered in the middle of the chain and zero at the 
boundaries. 

Recent applications of this method include the calculation of excitations in spin-orbital models 
(SU(4)) in a magnetic ficld[57], spin dynamics in models for cuprate spin ladders including cyclic 
spin exchange[75], optical conductivity of the ionic Hubbard model[248], excitations in the one- 
dimensional Bosc- Hubbard model[249] and the optical response in ID Mott insulators [250]. 

In this section we have presented a method of calculating dynamical responses with DMRG. 
Although the basis truncation is big, this method keeps only the most relevant states and, for 
example, even by considering a 0.1% of the total Hilbert space (for N = 28 only ~ 40000 states 
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Figure 4: Momentum weights of a target state with q = it for N = 28 (circles), AT = 44 (squares), 
N = 60 (diamonds) and N — 72 (triangles). The dotted lines are a guide to the eye. 

are kept) a reasonable description of the low energy excitations is obtained. We show that it is also 
possible to obtain states with well defined momenta if the appropriate target states are used. 

5.1.1 Correction vector technique 

Introduced in Ref. [251] in the DMRG context and improved in Ref. [247], this method focuses on 
a particular energy or energy window, allowing a more precise description in that range and the 
possibility of calculating spectra for higher energies. Instead of using the tridiagonalization of the 
Hamiltonian, but in a similar spirit regarding the important target states to be kept, the spectrum 
can be calculated for a given z = w + it] by using a correction vector (related to the operator A that 
can depend on momentum q). 

Following (6), the (complex) correction vector \x(z)) can be defined as: 

\x(z)) = j^A\^) (13) 

so the Green's function can be calculated as 

G{z) = {M^\x{z)) (14) 
Separating the correction vector in real and imaginary parts \x(z)) — \x r (z)) + i\x l (z)) we obtain 

{(H-wf+ V 2 )\x\ Z )) = -r ] A\^) (15) 

and 

\x r {z))= 1 -{w-H)\x i {z)) (16) 

The former equation is solved using the conjugate gradient method. In order to keep the infor- 
mation of the excitations at this particular energy the following states are targeted in the DMRG 
iterations: The ground state the first Lanczos vector A\i[)o) an d the correction vector \x{z)). 
Even though only a certain energy is focused on, DMRG gives the correct excitations for an energy 
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range surrounding this particular point so that by running several times for nearby frequencies, an 
approximate spectrum can be obtained for a wider region [247] . 

A variational formulation of the correction vector technique which leads to more accurate excited 
energies and spectral weights has been developed in [252] . It has been successfully applied to calculate 
the optical conductivity of Mott insulators [253]. 

5.2 Moment expansion 

This method [254] relies on a moment expansion of the dynamical correlations using sum rules that 
depend only on static correlation functions which can be calculated with DMRG. With these mo- 
ments, the Green's functions can be calculated using the maximum entropy method. 

The first step is the calculation of sum rules. As an example, and following [254], the spin-spin 
correlation function S z (q, w) of the Hciscnberg model is calculated where the operator A of Eq. (4) 
is S z (q) = N~ x / 2 S z (l) exp(iql) and the sum rules are[255]: 

f°° dwS z {q 7 w) 1 , 

J a IT w 2 

f°° dw S z (q,w) 1 . 
m 2 (q) = / — w w ' ' = -S'(q,t = 0) 
Jo " w 2 

mM = r-^^^ = - 1 -([[H,S z (q)],S z (- q )]) 
io " w 2 

= 2[l-co S (q)]J2(S+Sr +1 +Sr S + +1 ) (17) 

i 

where x(<Z, w — 0) is the static susceptibility. These sum rules can be easily generalized to higher 
moments: 

r°°dw l _ l SHvM 

Jo n w 
= -±([[H,...,[H,S'(q)U,S'(-q)]) (18) 

for I odd. A similar expression is obtained for I even, where the outer square bracket is replaced by 
an anticommutator and the total sign is changed. Here H appears in the commutator I — 2 times. 

Apart from the first moment which is given by the static susceptibility, all the other moments 
can be expressed as equal time correlations (using a symbolic manipulator) . The static susceptibility 
X is calculated by applying a small field h q rii cos(qi) and calculating the density response (n q ) = 
l/-^X)i( n «) c °s(qi) with DMRG. Then \ = ( n q)/h q for h q — > 0. These moments are calculated for 
several chain lengths and extrapolated to the infinite system. Once the moments are calculated, the 
final spectra is constructed via the Maximum Entropy method (ME), which has become a standard 
way to extract maximum information from incomplete data (for details see Ref. [254] and references 
therein). Reasonable spectra are obtained for the XY and isotropic models, although information 
about the exact position of the gaps has to be included. Otherwise, the spectra are only qualitatively 
correct. This method requires the calculation of a large amount of moments in order to get good 
results: The more information given to the ME equations, the better the result. 



5.3 Finite temperature dynamics 

In order to include temperature in the calculation of dynamical quantities, the Transfer Matrix 
RG described above (TMRG[185, 187, 189]) was extended to obtain imaginary time correlation 
functions[256, 257, 258]. After Fourier transformation in the imaginary time axis, analytic contin- 
uation from imaginary to real frequencies is done using maximum entropy (ME). The combination 
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of the TMRG and ME is free from statistical errors and the negative sign problem of Monte Carlo 
methods. Since we are dealing with the transfer matrix, the thermodynamic limit can be discussed 
directly without extrapolations. However, in the present scheme, only local quantities can be calcu- 
lated. 

A systematic investigation of local spectral functions is done in Ref. [258] for the anisotropic 
Heisenberg antiferromagnetic chain. The authors obtain good qualitative results especially for high 
temperatures but a quantitative description of peaks and gaps are beyond the method, due to the 
severe intrinsic limitation of the analytic continuation. This method was also applied with great 
success to the ID Kondo insulator [257]. The temperature dependence of the local density of states 
and local dynamic spin and charge correlation functions were calculated. 

6 Conclusions 

We have presented here a very brief description of the Density Matrix Rcnormalization Group 
technique, its applications and extensions. The aim of this article is to give the unexperienced 
reader an idea of the possibilities and scope of this powerful, though relatively simple method. The 
experienced reader can find here an extensive (however incomplete) list of references covering most 
applications to date using DMRG in a great variety of fields such as Condensed Matter, Statistical 
Mechanics and High Energy Physics. 
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